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,We study ultracold fermionic atoms trapped in a three dimensional optical lattice by combining the real-space dynamical mean-field 
approach with continuous-time quantum Monte Carlo simulations. For a spin-unpolarized system we show results the density and 
pair potential profile in the trap for a range of temperatures. We discuss how a polarized superfluid state is spatially realized in the 
spin-polarized system with harmonic confinement at low temperatures and present the local particle density, local magnetization, 
and pair potential. 
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1. Introduction 

Since the successful realization of Bose-Einstein condensa- 
tion in bosonic **^Rb and ^^Na J^lsystems, ultracold atoms 
have attracted considerable interest [Sl |4|, |5|] . One of the most 
active research areas in this field concerns optical lattice sys- 
tems, which are generated by subjecting the trapped ultracold 
atoms to aperiodic potential generated by appropriated laser 
beams H iSJ. The setups provide clean quantum sys- 
tems with parameters which can be tuned in a controlled fash- 
ion. Remarkable phenomena have been observed such as the 
phase transition between a Mott insulator and a superfluid in 
bosonic systems ifioll . In the fermionic case, both the super- 
fluid state ifTlll and the Mott insulating state 1 12, 13] have been 
observed. Furthermore, spin imbalanced populations have re- 
cently been realized 11511 . which stimulate further theoreti- 
cal and experimental investigations on ultracold fermionic sys- 
tems, and allow to investigate well-known ideas from conven- 
tional condensed matter physics. 

, For the imbalanced system with attractive interactions, inter- 
esting ordered ground states have been proposed as the stan- 
dard s-wave pairing at the Fermi surface is then modified. One 
of the most prominent candidates is the Fulde-Ferrell-Larkin- 
.Ovchinnikov (FFLO) phase |3 17], in which Cooper pairs 
■with nonzero total momentum are formed. This phase has been 
observed in the high field region in CsCoIns lialmM, and 
has theoretically been discussed in this compounds |2lll. as well 
as cold atoms with imbalanced populations 1221 12311 . Another 
proposed phase is the breached-pair (BP) phase, where both 
the superfluid order parameter and the magnetization are finite 



at zero temperature 1124 125L l26l 1271 12811 . When one considers 



three dimensional optical lattice systems at finite temperatures, 
the naively expected polarized superfluid state may be more sta- 
ble than the others. However, it is not clear how the polarized 
superfluid state is realized, and how the pair potential and the 



magnetization are spatially distributed in the imbalanced sys- 
tem with a confining potential. Understanding this issue may 
be important to observe the polarized superfluid state experi- 
mentally. 

In order to clarify these aspects, we investigate the attractive 
Hubbard model with imbalanced spin populations and a confin- 
ing potential. This allows to discuss the effect of the imbalanced 
spin populations on the superfluid state . By using real-space dy- 
namical mean-field theory (R-DMFT) ll29t[30i[3ll[32ll. we study 



the low temperature properties of the model. Here, we use the 
continuous-time quantum Monte Carlo (CTQMC) method ll33ll 
based on the Nambu formalism as an impurity solver ll34ll . By 
calculating the local particle density, local magnetization and 
pair potential, we clarify how the polarized superfluid state is 
realized in the spin imbalanced system within the confining po- 
tential. 

The paper is organized as follows. In Sec. |2] we introduce 
the model Hamiltonian and summarize various aspects of the 
R-DMFT. We demonstrate how the superfluid state is realized 
in a fermionic optical lattice with the confining potential in Sec. 
[3] A brief summary is given in the last section. 



2. Model Hamiltonian and Method 

Let us consider ultracold fermionic atoms in an optical lattice 
with harmonic confinement, which may be described by the fol- 
lowing attractive Hubbard Hamiltonian ll3l[36i[37i[38i[39il40[ 



<ij)<r 



(1) 
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where c,o-(c.^) annihilates (creates) a fermion at the /th site with 
spin 0" and - c\^Cja-. h acts as a magnetic field which allows 
as to tune the spin population imbalance, is the chemical po- 
tential, f(> 0) denotes the nearest neighbor hopping and U{> 0) 
the attractive interaction. V is the curvature of the harmonic po- 
tential. The notation {ij) indicates that the sum is restricted to 
nearest neighbors, r, is the distance measured from the center 
of the trap and a is the lattice spacing. 

The ground-state properties of the Hubbard model on in- 
homogeneous lattices have been studied theoretically by vari- 
ous methods such as the Bogoljubov-de Gennes (BdG) equa- 
tions 1I43I1 . the Gutzwiller approximation |44], the slave-boson 
mean-field approach 14511 . variational Monte Carlo simula- 
tions ll46ll . and the local density approximation! 47 1. On the 
other hand, there are few studies addressing the effect of imbal- 
anced populations beyond the static mean-field approach. The 
density matrix renormalizatio n group method 1481 149|] and the 



5111 are powerful for low di- 



quantum Monte Carlo method JSC 
mensional systems, but it may encounter difliculties when ap- 
plied it to higher dimensional systems. Here we use the R- 
DMFT approach il l32|], where local particle corre- 

lations are taken into account precisely. This treatment is for- 
mally exact for the homogeneous lattice model in infinite di- 
mensions H Sim in and the method has successfully been 
applied to some inhomogeneous correlated systems such as the 
surface lf52ll or the interface of a Mott insulators |53], and to 
fermionic atoms 1 54 , 55 , 5^ . 

In R-DMFT, the lattice model is mapped to a collection of ef- 
fective impurity models. The lattice Green function is then ob- 
tained via a self-consistency condition imposed on the impurity 
problems. When one describes the superfluid state in the frame- 
work of R-DMFT ililMIIillll, the lattice Green's function 
should be represented in the Nambu-Gor'kov formalism. For a 
system with L lattice sites it is then given by an (L x L) matrix, 
where each component consists of a (2 x 2) matrix as, 

%(icj„)] , (2) 

where = 1, 2, ■ • ■ , L, yu, = yL(-y(r,7a)^, iXj. is the z component 
of the Pauli matrix, &o the identity matrix, w„ - (2« + 1 )7tT the 
Matsubara frequency, and T the temperature, d^ij) is 1 when 
site / and j are neighboring sites and zero otherwise. The site- 
diagonal self-energy at the /th site is given by the following 
(2 X 2) matrix. 



2, {ico„) = 



£,|(/W„) SiiioJn) 

Si(ia)„) -2^(;w„) 



(3) 



where 'La-{ia>„) [S (ia>n)] is the normal (anomalous) part of the 
self-energy. In R-DMFT, the self-consistency condition is given 
by 



(4) 



rimp,i {i(^n) , 

where G„„p,,- is the Green's function of the eff'ective impurity 
model for the ;th site. Then the effective medium for each site 
is given by 

1-1 



In this paper, we focus on the low energy state without lattice 
symmetry breaking ||56| . In this case, the point group symmetry 
can be employed to efficiently deduce the lattice Green's func- 
tion. The inverse lattice Green's function eq. (|2]i is transformed 
in terms of a unitary matrix U, as 










Ma,, 














Mt2„ ) 



(6) 



where M, is the (m, x nii) matrix, »z,(< L) is the number of 
elements, and / runs over the representations of the group Oh- 
The lattice Green's function is then obtained by 



k m,n 


















[Mt,J 



U(kmU [Mi.^]m„ U(kn),i 



u, 

(7) 
(8) 



When one considers a system size r, < la, the total num- 
ber of sites L - 1419. In this case, there are 58 inequivalent 
sites. Therefore, by solving fifty-eight kinds of effective impu- 
rity models and making use of eq. ([8]l, we can iteratively solve 
the R-DMFT equations and discuss the stability of the polarized 
superfluid state in the large cluster. 

When R-DMFT is applied to our inhomogeneous system, 
it is necessary to solve a large number of effective impurity 
models. There are various numerical techniques such as exact 
diagonalizationjSTi and the numerical renormalization group 
[42, 58, 59, 60|. One of the most powerful methods is CTQMC, 
which has been developed recently. In this method, Monte 
Carlo samplings of collections of diagrams for the partition 
function are performed in continuous time. Therefore, the Trot- 
ter error, which originates from the Suzuki-Trotter decompo- 
sition, is avoided. Furthermore, this method is efficient and 
applicable to more general classes of models than, for exam- 
ple, the Hirsch-Fye algorithm l6lll . The CTQMC method has 
successfully been a ppli ed to various systems such as the Hub- 
bard model i34l62l 16311. the periodic Anderson model 116411 . the 
Kondo lattice model lisll . and the Holstein-Hubbard model. ll66ll 
Here, we use the continuous-time auxiUary field (CTAUX) ver- 
sion of weak-coupling CTQMC 134 16711 to discuss the super- 
fluid state in the optical lattice system. Some details of the 



CTQMC method are explained in Appendix A 



0i ' (ioj,,) 



i (ioj,,)] + % {ia>„) . 



(5) 



3. Results 

We now consider a three-dimensional optical lattice system 
with a confining potential to discuss how the superfluid state is 
realized in this spatially inhomogeneous set-up. In this paper, 



we use the hopping integral r = 1 as the unit of the energy 
and fix the confining potential as V = 0. 1 and the total particle 
number as N = 80. We calculate site-dependent static physical 
quantities such as the particle density «,-, the pair potential A,- 
and the magnetization m,-, which are defined by 



A,- = (c,TC,i) = F,(0+), 



(9) 

(10) 
(11) 



By performing R-DMFT with the CTQMC method, we obtain 
these quantities first for the balanced system {Ni = Ni - 40), 
as shown in Fig. [1] At high temperatures, particle correlations 
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Figure 1: Profiles of particle density («,o-) and pair potential A; as a function 
of r when Vq = OA, U = 8 and = Ni = 40. 



are small and kinetic energies high such that the fermions are 
widely distributed in the trap, as shown in Fig. [T](a). Decreas- 
ing the temperature, the particles gather around the center of 
the system, which is due to the existence of the attractive in- 
teraction. Below a certain critical temperature Tc - 0.4 ~ 0.5, 
the pair potential A, becomes finite, as shown in Fig. [1] (b). 
This means that a superfluid state is reaUzed in the region with 
Hi + 0. It is also found that the magnitude of pair potential 
depends on the site. This originates from the existence of the 
harmonic confinement in the system, which is consistent with 
the results obtained from BdG equations 143 1. 



We next discuss the effect of the imbalanced populations. 
The obtained results for P - 0.25 are shown in Fig. |2] where 
P - {Ni - Ni)/(Ni + Ni)is the spin imbalance parameter. This 
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Figure 2: Profiles of pair potential A, and local magnetization »!, as a function 
of r when Vq = 0. 1 , {7 = 8, A?^ = 50 and Ni = 30. Inset of (a) shows the profile 
of particle density (iij-f). 



case is similar to a system with applied magnetic field. There- 
fore, the superfluid state should become unstable 13411 . In fact, 
in contrast to the balanced system, the normal metallic state is 
stabilized even atT - 0.33. Further decrease of the temperature 
(r = 0. 1) induces the pair potential in a certain region (r,- < 4a), 
as shown in Fig. |2](a). Since the magnetization is finite in this 
region, we conclude that the polarized superfluid state is real- 
ized in the region (0 < r, < 4a). On the other hand, the normal 
metallic state still remains in the region (4a < r,- < 6a). Note 
that at r - 0.1, the pair potential has a maximum in the cen- 
ter of the system, where particle pairs are strongly formed. As 
can be seen in Fig. |2](b), the magnetization at the center of the 
trap first increases on reducing the temperature. However, when 
the system becomes superfluid, the magnetization is suppressed 
and the maximum is pushed to larger values of r,- resulting in a 
non-monotonic behavior appears in the magnetization curve. It 
is expected that further decrease of the temperature will lead to 
a complete suppression of the magnetization near the center of 
the trap. Detailed calculations will be presented elsewhere. 
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4. Summary 

We have investigated ultracold fermionic atoms trapped in an 
optical lattice with spin imbalanced populations. By combining 
the real-space dynamical mean-field theory with continuous- 
time quantum Monte Carlo simulations based on the Nambu 
formalism, we have calculated the local particle density, local 
magnetization, and the pair potential in the system. We have 
demonstrated how a polarized superfluid state is spatially real- 
ized at low temperatures in the model with harmonic confine- 
ment. 

In this paper, we have studied the stability of the polarized 
superfluid state in a system with small number of fermions. If 
the total particle number is larger, the local particle density may 
increase beyond «, - 0.5. It is known that the density wave state 
and the superfluid state are degenerate in the homogeneous sys- 
tem at half filling (« = 0.5) except in one dimension, which 
means that a supersolid state might be realizable in an optical 
lattice system with modulated particle density. It is an interest- 
ing problem to clarify this by means of our method. Work along 
those lines is in progress. 
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Appendix A. Continuous-Time Quantum Monte Carlo 
simulations in the Nambu Formalism 

This appendix explains some details of the CTAUX method 
1 6711 . When the superfluid state is discussed in the framework of 
DMFT, the total particle number is not conserved in the effec- 
tive Anderson impurity model. The Hamiltonian is then given 
by 



H - H() + Hjj, 

i/o = ^ epo-npa + ^ {Vpo-dla,,^ + li.c.^ 

per per 



Hu = -U 



ndf^di - ^{'U^+'idi - 1) 



(A.l) 

(A.2) 
(A.3) 



where a^o- (do-) annihilates a fermion with spin cr in the pth 
orbital of the effective baths (the impurity site). The effective 
bath is represented by e^o- and Ap, and Vpa- represents the hy- 
bridization between the effective bath and the impurity site. Eda 
is the energy level for the impurity site, n^o- - a^p^^apa-, and 



ndcr - do-da-- The Green's function should be defined by the 
2x2 matrix, as 



G(t) = 
where 



Gt(t) FiT) 
F*(t) -Gii-T) 



gat) = <r,c^(T)c;(0)>, 

F(t) = {TrC^(T)cim, 
F*(t) = {Trcl(T)c\{0)), 



(A.4) 



(A.5) 
(A.6) 
(A.7) 



where Tr is the imaginary-time ordering operator and we have 
chosen the Green's functions Go-(t) to be positive. 

To perform simulations, we consider here a weak coupling 
CTQMC approach. The partition function Z is given by 



Z = Tr e 



Tre 



-jl^dTH.jT) 



n=0 ^" •Jt„-\ 

X (-l)«Tr[e-^«'//2(T„)//2(T„-i)---//2(Ti)], (A.8) 

where H2(j) - e^^^H2e^'^^' and yS = l/T. Here, we have di- 
vided the impurity Hamiltonian Eq. ( lA.lb into two parts as. 



H\ — H - Hi, 

H2 = Hu-Kip 
K 



= - y 



,yj(n,+;H-l) 



(A.9) 



(A. 10) 



with y - cosh '(1 + /3U/2K), and K some nonzero constant. 
In this paper, we set A' = 1 in the CTQMC simulations. The 
introduction of the Ising variable s in H2 enables us to per- 
form simulations at arbitrary filling. An «th order configura- 
tion c - {ii, S2, ■Sn; Ti, 72, T„) corresponding to auxil- 
iary spins si,S2, - ■ - , s„ at imaginary times ri < T2 < . . . < t„ 
contributes a weight 



(^)^-''2''Zodet[^(«)]-' (A.11) 



to the partition function. Here, Zq - Tr[e ] and .^V*"' is an 
n X « matrix, where each element consists of a 2 x 2 matrix: 



' ^ f («) _ |(«) (n) _ /(n)^ ^ 



In' = Sijo-Q, 



ij 

pin) 
•J 



s, 



(«) 



gO]{Ti - Tj) foiTi - Tj) 
-fo('^i-'^j) 80l(Tj-Ti) /' 



(A. 12) 
(A.13) 
(A. 14) 

(A. 15) 



with i,j = 1,2, •■■«. The sampling process must satisfy er- 
godicity and (as a sufficient condition) detailed balance. For 
ergodicity, it is enough to insert or remove the Ising variables 
with random orientations at random times to generate all possi- 
ble configurations. 
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To satisfy the detailed balance condition, we decompose the 
transition probability as 



P (i j) 



(A. 16) 



where pP"P(p'"^'^) is the probability to propose (accept) the tran- 
sition from the configuration ; to the configuration j. Here, we 
consider the insertion and removal of the Ising spins as one step 
of the simulation process, which corresponds to a change of + 1 
in the perturbation order The probability of insertion/removal 
of an Ising spin is then given by 



^Pl"op 



(n 



dT 

n + \) - — , 
' 2/3 



'(« + 1 ^ «) = 



1 



(A. 17) 



(A.18) 



« + 1 

For this choice, the ratio of the acceptance probabilities be- 
comes 

, detN^"'> 



■in 



1) 



K 



pace („ + 1 ^ „ + 1 



detM"+i'' 



(A. 19) 



When the Metropolis algorithm is used to sample the configu- 
rations, we accept the transition from n to « + 1 with the proba- 
bility 



imn 



1, 



pace + 1 _> 



(A.20) 



In each Monte Carlo step, we measure the following Green's 
functions (0 < r < /?), 



GAr) 
F(r) 
F*(t) 



Z 

Z L 

= iTrf, 
7 L 



Z L I 



CT(T)Ci(0)l , 



(A.21) 
(A.22) 
(A.23) 



By using Wick's theorem, the contribution of a certain configu- 
ration c is given by 



G^(t) = det[A?(")]det 
F'(t) = det[A?("']det 
F*'iT) = det[M"']det 













Q' 




Mt) 




G" 


^ R*' 





(A.24) 
(A.25) 
(A.26) 



where Qa-, Q', Q*',Ra,R',R*' are vectors, in which the /th ele- 
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